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The emergence of the Haldane Chern insulator state due to strong short range repulsive interac¬ 
tions in the half-filled fermionic spinless honeycomb lattice model has been proposed and challenged 
with different methods and yet it still remains controversial. In this work we revisit the problem 
using the infinite density matrix renormalization group method and report numerical evidence sup¬ 
porting i) the absence of the Chern insulator state, ii) two previously unnoticed charge ordered 
phases and iii) the existence and stability of all the non-topological competing orders that were 
found previously within mean held. In addition, we discuss the nature of the corresponding phase 
transitions based on our numerical data. Our work establishes the phase diagram of the half-filled 
honeycomb lattice model tilting the balance towards the absence of a Chern insulator phase for this 
model. 


I. INTRODUCTION 


Topological phases of matter are remarkably robust 
states; their responses to external fields are governed by 
topological invariants and thus many of their most impor¬ 
tant properties are insensitive to local perturbations. 2 3 
Topologically protected responses result in intrinsically 
novel phenomena including fractionalization of quantum 
number^ or transport governed by quantum anomalies 5 
and can lead to diverse applications, ranging from fault 
tolerant quantum computation to spintronicsP® Added 
to the remarkable experimental discoveries of new topo¬ 
logical phases in both two and three dimensions, these 
ideas have boosted a sustained and voluminous scientific 
effort in the last decade that attempts to classify them 
and determine when and how can they emerge P 
This work addresses a concrete question regarding the 
emergence of topological phases that has so far remained 
controversial: the existence of the Chern insulator phase 
triggered by short range repulsive interactions in the 
fermionic spinless half-filled honeycomb tight binding 
model. The Chern insulator state, first identified by 
Haldane^ in the particular case of the honeycomb lattice, 
is a zero field analogue of the integer quantum Hall 
effect; the Hall conductivity contribution of each band 
is quantized in integer units of e 2 /h and determined by 
a topological invariant, the Chern number. Realising a 
Chern insulator in nature is not a trivial task; 8 time- 
reversal must be broken with the resulting magnetic flux 
averaging to zero over the unit cell, and hence over the 
entire sample. 

In a proof of principle, a set of mean field studied®^ 
discussed how the Chern insulator can emerge from 
repulsive short range interactions. It was shown that 
spinless fermions hopping in the half-filled honeycomb 
lattice with nearest and next-to-nearest neighbor inter¬ 
actions, Vi and V 2 respectively [see Fig. [I], displayed a 
Chern insulator phase as described by Haldane P The 


non-interacting realization of this Chern insulator phase 
has complex next-to-nearest neighbor hoppings [see 
Fig. [2jg)]. Within the mean field paradigm it occurs in a 
region with V 2 > Vi where V 2 spontaneously breaks time 
reversal symmetry generating complex next-to-nearest 
neighbor hopping strengths. The latter condition was 
shown not to be a generic requirement; Chern insulator 
phases can be realized in mean field with only V± at 
the ex pense of enlarging the unit cell and doping the 
system. 12 13 In addition, analogous topological phases 
have been obtained in other lattices such as the 7r— flux 
modethHASJ an d the Kagome latticed 

The important question that still remains to be an¬ 
swered is whether the interaction induced Chern state 
survives after the effect of quantum fluctuations is in¬ 
cluded. The Haldane Chern insulator phase competes 
with more conventional but also interesting orders, that 
can jeopardize its emergence and are depicted schemat¬ 
ically in Fig. [2] In the V\ V 2 limit, a charge ordered 
state depicted in Fig. [2jb) is expected, where the A and 
B sublattices are populated differently P Being con¬ 
nected to the classical ground state at V\/t —> oo with 
only one occupied sublattice, this state (CDW I) has 
been found to be very robust both in mean field and 
exact diagonalization. 18 1 -® For V 2 V± on the other 
hand it was shown under the mean field paradigni 12 ^ 
that the phase space originally attributed to the Hal¬ 
dane Chern insulator was severely reduced by the pres¬ 
ence of a sublattice charge modulated state (or CMs) 
depicted in Fig. [2jc) . The CMs phase was corrobo¬ 
rated to survive quantum fluctuations in exact diag¬ 
onalization 18 ~ 20 and within a variational Monte Carlo 
approach. 2 ® The exact nature of the phase was further 
discussed in Ref. QJ5] by computing the charge structure 
factor; its suppressed Fourier component at T suggested 
the absence of charge imbalance between the sublattices 
and was therefore termed CM phase. In addition, at in¬ 
termediate Vi > V 2 a Kekule bond ordei^HH] depicted 
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diagonalization,^ further insights are needed to cor¬ 
roborate its existence in the thermodynamic limit. 
Finally, the Haldane Chern insulator phase of Fig. 0(g), 
stable within mean field, was found however to be 
absent in exact diagonalization with periodic boundary 
conditional 19 and cluster perturbation theorjDH sug¬ 
gesting that quantum fluctuations indeed can destabilize 
its emergence. Nonetheless, this interpretation was 
challenged by Ref. [20] that observed hints of this phase 
in exact diagonalization with open boundary conditions 
and variational Monte Carlo. 


FIG. 1. The top left illustration shows a schematic represen¬ 
tation of the interacting tight binding model considered in 
this work: fermions hopping in a half-filled honeycomb lattice 
with nearest-neighbor hopping t (solid lines) and interacting 
via nearest- and next-to-nearest neighbor interactions, Vi and 
Vi respectively (curved lines) as defined by 0 - The right il¬ 
lustration shows the iDMRG unit cell used in this work with 
3x6 unit cells yielding a cylinder with a circumference of 
L y = 12 sites depicted in the bottom left illustration. 


in Fig.[2jd) was found to be stable within mean field 
theory!® 1 ' 

This three-fold degenerate bond order triples the 
original honeycomb two atom unit cell breaking its 
translational symmetry. Although numerical evidence 
consistent with such bond order was found in exact 



(e) (f) (g) 


FIG. 2. Different orders allowed in the half-filled spinless hon¬ 
eycomb lattice considered in this work: (a) semimetal, (b) 
charge density wave I (CDW I), (c) sublattice charge mod¬ 
ulated (CMs), (d) Kekule, (e) CDW II, (f) CDW III and 
(g) Haldane Chern insulator phases. Phases (a)-(d) and (g) 
were found within mean field theory by Refs. [9] HOI and 1T21 
The survival of the Haldane Chern insulator phase (g) was 
challenged within exact diagonalization with periodic bound¬ 
ary conditions by Refs. [T 8 | and [19] but found in Ref. [20] for 
open boundary conditions. The generic charge imbalance be¬ 
tween the sublattices in phase (c) was not found within exact 
diagonalization . 19 The two sublattices are depicted in blue 
and orange, 0 < A < 0 < 1/2 with g + A < 1/2 describe the 
deviations from half-filling per site, i.e. (n) = 1/2 ± a with 
a € {g,g± A}. 


The contradictory numerical evidence regarding the 
existence of the Chern insulator phase in particular, 
and other competing phases in general, needs of an 
alternative approach that includes quantum fluctua¬ 
tions while minimizing finite size effects. In this work 
we therefore revisit the controversies left unsolved 
by previous studies using the infinite density matrix 
renormalization group method (iDMRG/^^II. This 
variational method determines the ground state of sys¬ 
tems of size L x x L y where L x is in the thermodynamic 
limit and L y goes beyond what is achievable in exact 
diagonalization. Traditionally a method for finding the 
ground state of one-dimensional systems, iDMRG has 
recently been successfully applied to two-dimensional 
systems. The infinite and finite DMRG method were 
indeed shown to be able to characterize the properties 
of fractional quantum HallpSl Z 2 quantum spin liquid,^ 
chiral spin liquid 30 and bosonic and fermionic fractional 
Chern insulating states.® 1 -®^ As long as entanglement 
remains low and the state has short correlation lengths 
the ground state can be represented faithfully by a 
product of matrices -termed matrix product state 
(MPS)- of a computationally affordable dimension y, 
known as the bond dimension PSHUI This requirement 
is commonly met by gapped systems and in particular 
by the orders argued above to be expected instabilities 
in the half-filled honeycomb lattice of interest here. 
Moreover, the iDMRG method deals with an infinite 
system and thus, given a suitable size of the unit cell 
used for the simulations, it can directly probe ground 
states that spontaneously break the original symmetries 
of the hamiltonian and the phase transitions among 
them, while allowing for quantum fluctuations to play a 
role. 

Motivated by these advantages we have mapped out 
the {V[, V 2 } phase diagram for the half-filled honeycomb 
lattice in the infinite cylinder geometry with the iDMRG 
method [see Fig. 0. Our main findings are summarized 
next. First we show compelling numerical evidence that 
supports the absence of the Chern insulator state in a 
wide region of { V \, V 2 } phase space. Second, we find two 
new p hase s not reported previously neither within mean 
fielcP 10 12 nor the first exact diagonalization resultd^® 
and analyze their semiclassical features. Third, we 
provide numerical support for the existence and stability 
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of all the competing orders that were found previously 
within mean field. In particular we characterize the 
two more controversial states, the CMs and Kekule 
states, by computing bond and charge expectation 
values and present further arguments of the existence 
of the former from the semiclassical large interaction 
limit. Fourth, we provide numerical evidence to assess 
the first or second order character of the relevant phase 
transitions. Concretely we address this issue by identi¬ 
fying the features that the correlation length £ and the 
entanglement entropy S present as a function of {V}, 14}- 

This work is structured as follows. After describing 
the method and model in section m we characterize the 
different phases that we find in section |III| The corre¬ 
sponding phase transitions among them are analyzed in 
section |IV| and we end with a discussion and prospect of 
our results in section [V] Appendix [A] includes details of 
our semiclassical analysis as well a discussion of the CMs 
state order parameter. 

II. MODEL AND METHOD 

We investigate a system of spinless fermions hopping 
on a honeycomb lattice with real nearest neighbor hop¬ 
ping t > 0 interacting via nearest and next-to-nearest 
neighbor interactions {Vi ,V '2 } > {0,0} respectively [see 
Fig. 0 - The Hamiltonian for this system can be written 
as 

H = —t y~](cjc,-+h.c.)+Vi y: njUj + V 2 ^ riiUj, (1) 

(*d> (»d> <(hi» 

where c i (c|) annihilates (creates) an electron at the Ath 
site of the honeycomb lattice. 

In order to find the ground state of the system in 
the { V [. V 2 } phase space we employ the iDMRG algo¬ 
rithm —" 27 on an infinite cylinder geometry [see Fig. Fy. 
As discussed above, such infinite geometry allows to 
probe ground states with spontaneous symmetry break¬ 
ing while taking into account quantum fluctuations. Fea¬ 
sible system sizes for our purposes are cylinders of cir¬ 
cumference L y = 6 , 8 ,10 and 12, given the exponential 
growth of computational cost with the circumference. 
For our calculations we choose L y = 12 and a unit cell of 
36 sites depicted in Fig. [T] 

We pick this geometry for the following two reasons. 
First, we have to keep in mind the structure of the recip¬ 
rocal space. Since we work in the thermodynamic limit 
in the ^-direction along the cylinder, the momentum in 
x-direction is a continuous quantity. However in the fi¬ 
nite y-direction around the cylinder, the momentum is a 
discrete variable. In the non-interacting case, the model 
forms a Dirac semimetal and the Fermi surface is located 
at the K and K' points.^ In order to capture the low 
energy physics correctly, it is of key importance that K 
and K' are allowed momenta in our unit cell which im¬ 
plies that only L y = 6 and 12 are suitable circumferences 


of our cylinder. Second, we have the freedom of fixing 
the length of our unit cell in the cc-direction since the 
computational for that only scales linearly. By choosing 
three rings of 12 sites each, all orders are commensurate 
with our unit cell and a further enlargement would not 
lead to different results in the parameter region we inves¬ 
tigated. All data we show in the remainder of the text 
is computed for L y = 12 and bond dimension \ = 1600 
unless otherwise stated. 


III. PHASE DIAGRAM 

We have mapped the phase diagram as a function of 
{ V [, V' 2 } with the method described above. Our main 
results are summarized in Fig. [3] We find six different 
phases: two new charge order phases labeled CDW II and 
CDW III and four previously reported mean field orders, 
the Kekule, CMs, CDW I and semimetal phases. 

We have characterized each phase through their charge 
and bond ground state expectation values. At each site 
the charge expectation value is defined by 

n* = (ctci). (2) 

The bond ground state expectation value on the other 
hand is defined as 

tij = (c\cj+ h.c.^. (3) 

Next we discuss the features of the different phases in the 
phase diagram in Fig. [3] 

A. Semimetal phase 

We start by discussing the semimetal phase in the 
phase diagram shown in Fig. [3j At Vi ft = Vi/t = 0 
the honeycomb lattice with nearest neighbor hopping its 
known to be described by a low energy theory in terms of 
two massless Dirac fermions! 35 Short range int eract ions 
are irrelevant in the renormalization group sensd 3 ® 37 and 
therefore they can only drive a transition to an ordered 
state when they have a magnitude comparable to the 
nearest neighbor hopping strength t. Such perturbative 
analysis guarantees that the semimetal is stable within 
the region {V\. 14} < t, only allowing for a uniform 
renormalization of the hopping strength t by interactions. 

For the semimetal phase in Fig. [3] and to numerical 
accuracy, we find by computing the charge expectation 
value ([ 2 ]) that m = 1/2 for all sites, indicating that this 
phase is not charge ordered. 

From © and choosing i,j to be nearest neighbors we 
find a small asymmetry between bonds pointing along 
and around the cylinder axis, with a relative difference of 
~ 10~ 2 . Such asymmetry should -up to a very small ef¬ 
fect due to the cylinder geometry which is always present 
for finite L y - vanish in a perfect semimetallic phase and 
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FIG. 3. Left: Phase diagram obtained with iDMRG calculations on an infinite cylinder of circumference L y = 12, \ = 1600. 
The phase boundaries, especially at second order transitions (see Sec. IV ) are to be taken with some error which can be 
estimated from Figs. [4] to [8] We draw sharp lines here for better clarity. Right: Representative charge and bond strength 
patterns for the four phases with largest unit cell discussed in the main text. The area of the blue circles is proportional to 
the particle number expectation value on the respective site given by Eq. |2|. The thickness of the ellipsoids on the bonds is 
proportional to the amplitude tij between nearest neighbors defined by Eq. ( 3J). The unit cells for each phase are depicted by 
the red polygons. These correspond to (a) CMs phase with Vi/t = 0.8, V 2 /t = 3.2, (b) Kekule phase with Vi/t = 5.6, Vi/t = 1.6, 
(c) CDW II phase with Vi/t = 5.6, Vi/t = 3.2, (d) CDW III phase with Vi/t = 9.2, Vi/t = 2.5. For (c) dashed lines indicate 
defect lines of rotated hexagons that have zero energetic cost in the classical limit (see main text). 


indeed it is severely reduced as the bond dimension \ is 
increased. This suggests that the cylinder geometry im¬ 
plemented in iDMRG artificially differentiates bonds in 
its two perpendicular directions. The asymmetry induced 
by the finite bond dimension vanishes as x is increased 
and thus it is not a physical effect. This is consistent 
with a semimetal state; the logarithmic divergence of en¬ 
tanglement of a metallic state requires a matrix product 
state with \ —► oo. The bond asymmetry reduces the 
entanglement by shifting the Dirac cones away from the 
K and K' points.®! 

Together, the previous numerical evidence are consistent 
with the semimetallic state. We note that numerically 
the semimetal state extends beyond {Vi, V 2 j < t towards 
higher interaction strengths through a narrow semimetal 
trench in the phase diagram. The larger size of this re¬ 
gion at L y = 6 (not shown) suggest that it is likely to 
shrink as the circumference of the cylinder is increased 
but a definitive statement requires going beyond the nu¬ 
merically accessible sizes. 

B. Charge density wave I (CDW I) 

Upon increasing Vi we identify a charge density wave 
state labeled CDW I in the phase diagram of Fig. [3| This 
state has a two site unit cell and is characterized by a 
finite order parameter of symmetry B 2 under the sym¬ 
metry group CgJ^ 

B 2 = (n A ) - ( n B ), (4) 

that can be calculated using |2]), where n A and n B are 
the fermion densities in the A and B sublattice sites of 
the two site unit cell respectively. The resulting charge 
order is depicted schematically in Fig. |2jb) . For instance, 


at V\/t = 4 and V 2 /t = 0.4 we find that B 2 = 0.817. The 
magnitude of B 2 increases with Vl and to numerical ac¬ 
curacy this phase has no appreciable bond order. 

For large Vi > i such a state is a natural insta¬ 
bility since the energy is minimized by a charge im¬ 
balance between the two sublattices. Inde ed, it has 
been found in a mean field approximation JUMUl ex¬ 
act diagonalizationp^Si an d quantum Monte Carlo 
simulations? 16 17 

C. Sublattice charge modulated phase (CMs) 

For small t -C V 2 the ground state is classically de¬ 
generate. Within mean field it was shown that as long 
as V 2 > Vi the system chooses a charge ordered pattern 
with charge modulation of wavevector K or K', termed 
the CMs phase®! and depicted in Fig. [2jc) . The order 
parameters for any such modulation can be taken as the 
corresponding Fourier components of the charge expecta¬ 
tion values of Eq. ([2]), namely (n^^K)) = (n AtB ( K'))*, 
where A and B indicate the two sublattices. It is con¬ 
venient to arrange these order parameters in a basis 
that has well defined transformation properties under the 
symmetry of the lattice. Under the action of the symme¬ 
try group Cq v , these order param eters transform as the 
four dimensional G' representation 39 ®! 


Gi a = ReMK)-n B (K)], 

(5) 

Gi B = Im[n A (K)-n B (K)], 

(6) 

G' 2x = lm[n A (K) + n B (K)], 

(7) 

G' 2y =Re[n A (K)+n B (K)]. 

(8) 


The CMs state corresponds to a finite expectation value 
of both the B 2 order parameter and G' = (1,0, 0, 0) and 
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the states generated from these by operations of the sym¬ 
metry group. A scalar order parameter for the CMs 
phase can be defined as S 2 Re[(G' 1 ) 3 — 3Gj (Gj) 2 ], where 
G[ = G' ix + iG' iy (see appendix A). The CMs structure 
is schematically shown in Fig. [2jc) . It minimizes the 
large cost attributed to V 2 by reversing two nearest neigh¬ 
bor dimers at the expense of paying the small energetic 
cost determined by Vj. Within exact diagonalization it 
has been argued that the CMs state survives quantum 
fluctuations BSE? 0 However, the sublattice charge imbal¬ 
ance predicted in mean field — was suggested to vanish in 
Ref. jTS] due to the absence of an enhanced Fourier com¬ 
ponent of the charge structure factor at the T point. We 
note that the CMs scalar order parameter defined above 
quantifies the charge imbalance since B 2 has to be non¬ 
zero for it not to vanish. 

With the iDMRG method we find that the CMs appears 
directly above the semimetallic phase [see Fig. [ 3 ]. A sam¬ 
ple of the six site unit cell charge and bond pattern ob¬ 
tained from Q for Vj jt = 0.8, V 2 /t = 3.2 is shown in 
Fig Ht a). The circles at each site have an area propor¬ 
tional to the strength of the charge at the given site, 
while the thickness of the links between the sites rep¬ 
resent the bond strengths. This bond and charge pat¬ 
tern survives in the region labeled CMs of the phase dia¬ 
gram and coincides with that obtained within mean field 
theorjEl depicted schematically in Fig. j2jc). The cor¬ 
responding charge values for Vj/i = 0.8, U 2 /t = 3.2 are 
given by g = 0.364 and A = 0.092. Moreover, we find a 
clear sublattice imbalance that grows as V 2 is increased 
that justifies the label CMs rather than the simpler CM. 
We note also that the bond ordering of this phase is as¬ 
sociated to the charge order: two neighboring sites with 
similar high (or low) charge densities that deviate from 
half-filling suppress the hopping between them due to the 
large (small) number of filled (empty) states at that site. 
We have also accounted for the existence of the CMs state 
from a strong coupling perturbation theory analysis. By 
exactly diagonalizing the interaction part of H in Eq. ([!]) 
and then finding the first order hopping corrections i/Vj 2 
in perturbation theory, we determined the ground state 
| GS) of a 3x 3 cluster. We then computed different scalar 
correlation functions of the charge order parameters, and 
found that only (GS'|B 2 Re[(G , 1 ) 3 - 3Gj(G' 2 ) 2 ]|GS) is fi¬ 
nite. This confirms that in the strong coupling limit, the 
ground state is indeed of the CMs form. The details of 
the procedure are presented in appendix [A] 


D. Kekule bond order 

The next phase we identify is the Kekule bond order. 
Like the CMs, it has a six site unit cell depicted schemati¬ 
cally in Fig. [2jd) with uniform charge order and two types 
of bonds, strong and weak. Under the mean field approx¬ 
imation the state arises between the CMs and the CDW I 
state. Although in previous exact diagonalization hints 
of this state are also observed, the evidence supporting 


its occurrence is not entirely conclusiv^^ and its tripled 
unit cell turns the analysis of larger clusters challenging. 

Within iDMRG we find that this state is stable in a 
finite but smaller region compared to both exact diago¬ 
nalization and mean field. In contrast to the anisotropy 
found in the semimetal phase, the Kekule order remains 
stable as the bond dimension \ is increased and stays 
finite upon extrapolation to infinite %. We show a rep¬ 
resentative of the state’s numerically obtained charge 
and bond order patterns in Fig. [ 3 ] (b) calculated with 
Eqs. ([2| and |3]) , respectively. The bond thickness is pro¬ 
portional to its strength. From Fig. [3](b) it is apparent 
that the state has indeed two bond strengths arranged 
as in Fig. §d) and no charge order. Therefore the sys¬ 
tem chooses to break the original two site unit cell trans¬ 
lational symmetry falling into one of the three distinct 
Kekule ground states.^ At Vj/t = 5.6 and V 2 /t = 1.6 the 
two types bonds are found numerically to be t\ = 0.522 
and t 2 = 0.428. 


E. Charge density wave II (CDW II) 

All phases that we have described so far do not differ 
from those predicted by mean field®E] anc [ Hie 
exact diagonalization studies.^’ 20 At finite Vj and 
sufficiently large V 2 we observe a different charge order 
state (CDW II), see Fig. [3] left panel. 

To gain a first insight on this phase before analyzing the 
numerical data it is instructive to discuss its classical 
limit. A strong coupling analysis (see appendix [A] for 
details) reveals quite generically that at t = 0 there are 
two classical regimes separated by the Vj = 4V 2 line. 
For Vj > 4U 2 the state CDW I is favored with the same 
charge pattern as in Fig. §b) but classical occupations 
0 or 1. At Vj = 4U 2 there exists a classically degenerate 
ground state manifold, to be discussed in the next sec¬ 
tion in the context of the CDW III phase. The analysis 
for Vj < 4U 2 with t = 0 reveals that there are essentially 
two types of degenerate states according to their trans¬ 
lational symmetry. The first is a stripe-like phase with 
a four atom unit cell depicted in Fig. [2je) (see also HIT) 
that is six fold degenerate in the thermodynamic limit. 
The second is a set of states that can be understood by 
taking the stripe-like phase and rotating 60° hexagons 
along a defect line that should wrap around the chosen 
cluster. Classically, the energy per site of these two types 
of states does not have first order hopping corrections 
and it is given by Ecdw ii = Vj/2 + Vj/4 + 0(t 2 /X j :2 ) 
(see appendix [A] for details). The first nontrivial order 
in perturbation theory depends on the chosen cluster; 
eventually the degeneracy between the six-fold degen¬ 
erate stripe-like phase and the phase with defects will 
be lifted, making one of them more favorable at finite 
The question then becomes which type of ground 
state does the system choose in the thermodynamic limit. 

Although we cannot attempt to fully answer this ques- 
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tion some light can be shed through our iDMRG numer¬ 
ical data. We indeed find that both configurations, with 
and without defects are ground states with very simi¬ 
lar energy, even close to the transition to the semimetal 
phase. A sample of one of this charge ordered patterns 
is shown in Fig. [3jc) for V\/t = 5.6 and Vjj/t = 3.2. This 
particular sample state has one defect line that wraps 
twice around the cluster (indicated by dashed lines); ro¬ 
tating the hexagons along these lines by 60° in an anti¬ 
clockwise direction we recover the more symmetric stripe- 
like configuration in Fig. [2je). As with the CMs state, we 
find that its bond order follows directly from its charge 
order, suppressed for neighboring sites that have equal 
filling. 

By initializing the algorithm with each inequivalent clas¬ 
sical ground state configuration it is possible to calculate 
to high precision each energy to then find, by compari¬ 
son, the lowest energy state at that particular 
We have determined that within the CDW II region the 
lowest energy is often achieved by a superposition of two 
classical ground states and their particle-hole conjugates. 
As expected from general arguments^ this state shows 
a degenerate low energy entanglement spectrum. Inspec¬ 
tion of the two point density-density correlation functions 
indicates that the more favored superposition of classical 
ground states are those which favor a uniform charge 
density pattern that at the same time distinguishes the 
bonds around and along the cylinder geometry. 


F. Charge density wave III (CDW III) 

At high Vj and between the previously discussed 
CDW II and CDW I states we find a third type of 
charge order with a twelve site unit cell labeled CDW III. 
A representative pattern of its charge and bond order 
as obtained numerically with © and ([3]) respectively 
is shown in Fig. §<*)■ As schematically represented in 
Fig- if), this state has three different occupations: half- 
hlling and 1/2 ± g. For the particular case of V\/t = 9.2 
and V' 2 /t = 2.5 we find that g = 0.223. Its bond or¬ 
der is determined by the charge occupations similar to 
the CMs and CDW II phases. As with the CDW II 
phase, th is phas e escaped identification in the original 
mean held® 10 12 and the first subsequent exact diagonal- 
izatiorP^SI studies due to its large unit cell. 

The emergence of the CDW III can be understood from 
a semiclassical point of view (see appendix A). At t = 0, 
the classical ground state of lowest energy for a cluster 
commensurate with this phase is either in the CDW II 
or CDW I class, and the two phases are separated by 
the line V\ = 414- Exactly at this line, however, a third 
classical state is degenerate in energy with these but, un¬ 
like the previous two, is affected by first order quantum 
corrections. This implies that around the line V\ ~ 4V2 
there is a finite strip of a new phase in the phase di¬ 
agram once finite t is included. We have checked that 
the two point correlation function calculated numerically 


from the iDMRG data and semiclassically agree qualita¬ 
tively, resulting in the charge distribution of the CDW III 
shown in Fig. [3](d). We find furthermore that these are 
consistent with those reported in the latest exact diago- 
nalization study.^ 


IV. PHASE TRANSITIONS 

In order to unravel the character of the corresponding 
phase transitions we now study two quantities that are 
sensitive to a phase change, the entanglement entropy S 
and the correlation length £. The entanglement entropy 
is defined by 

S = -Trp L lnp L , (9) 

where pl = Tr^|- j/>) ("01 is the reduced density matrix of 
the left semi-infinite half L of the cylinder after tracing 
out the right half R. The correlation length £ can be 
computed from the resulting state of the iDMRG cal¬ 
culation according to Ref. [57.. Both quantities are ex¬ 
pected to show a finite discontinuity when crossing a 
first order transition while the correlation length diverges 
with increasing bond dimension at a second order critical 
pointIn what follows, we focus on the lines labelled 
by cuts A-E in the phase diagram in Fig. [3j 

A. Cut A: semimetal—CDW I 

The first horizontal cut, labeled A addresses the char¬ 
acter of the phase transition between the semimetallic 
phase and the simplest charge density wave (CDW I) by 
fixing V 2 = 0. This phase transition has been pre viously 
addressed with the quantum Monte-Carlo methocP*® 7 ^ 
and was determined to be of second order character. The 
transition point with divergent correlation length was de¬ 
termined to be at Vi/£ = 1.356 via finite size scaling. In 
Fig-i we show the correlation length as a function of V\ 
for different values of the bond dimension 

Firstly, for V\ < 1.5t we observe that the correlation 
length drops as a function of Vj and diverges as the 
bond dimension \ is increased. This behaviour is ex¬ 
pected for a critical state such as the semimetal phase; 
the logarithmic divergence of entanglement of a metal¬ 
lic state requires an matrix product state with \ —» 00. 
For Vi > 1.5 1 the correlation length continues to drop 
but has no longer a significant dependence on This 
is characteristic of a gapped phase such as the charge 
density wave. In the inset of Fig. [4] we plot the entangle¬ 
ment entropy as a function of the logarithm of the bond 
dimension for two values of V) representative for the two 
phases. In the semimetal phase for VI = 0.8, we observe 
a finite entanglement scaling of S oc ln% characteristic 
of a critical phased whereas the entanglement entropy at 
V\ = 2.8 in the gapped CDW I phase is independent of 
the bond dimension. 
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FIG. 4. Correlation length at V 2 = 0, labeled cut A in Fig. [3] 
probing the transition between the semimetal and CDW I 
phases. The smooth behavior of £ indicates a second or¬ 
der transition. The dashed grey line at Vi/t = 1.356 shows 
where the transition has been detected in quantum Monte 
Carlo simulationsInset: Scaling of the entanglement en¬ 
tropy S with the bond dimension y f° r two values of Vi in 
the semimetal and CDW I phase. Unlike the CDW I phase, 
the semimetal phase presents the typical critical entanglement 
scaling S oc In y. 


The crossover between the two phases is smooth, sig¬ 
naling a second order phase tra nsitio n, in agreement with 
quantum Monte Carlo studies.® 17 With iDMRG it is 
however not possible to pin point the exact value of V\ 
where the transition happens via finite size scaling due 
to the few cylinder sizes we have available, as discussed 
in section [IT] However, from Fig. [4] it is possible to de¬ 
fine a crossover region of 1.3 < Vi 5. 1.5, consistent with 
the quantum Monte Carlo data,® 117 where the transition 
occurs. 


B. Cut B: CMs-CDW II 

In cut B we fix V 2 = 3.6 1 and study the phase transition 
between the CMs and the CDW II phases. The entangle¬ 
ment entropy is shown in Fig. [5] as a function of V\. This 
quantity shows a clear discontinuity at V 2 /t ~ 2.2 sig¬ 
naling a direct first order phase transition between these 
two phases. We observe the same behavior when looking 
at the correlation length not presented here. Addition¬ 
ally, the energy of the ground state as a function of V 2 
shows a kink at the same critical V 2 which further sup¬ 
ports the presence of a level crossing at that point. As 
expected for gapped phases, the entanglement entropy 
only depends very weakly on the bond dimension, in the 
CDW II phase even a very low y °f 400 is sufficient to 
faithfully represent the state. 

A first order phase transition is also expected from the 
strong coupling approach. As detailed in appendix [X] 
the energy per site of the CM phase in a 3 x 3 cluster is 


FIG. 5. Entanglement entropy at V 2 /t = 3.6, labeled cut B 
in Fig. [3] The finite discontinuity of S at V 2 /t « 2.2 clearly 
signals a first oder transition between the CMs and CDW II 
phases. 


Ecms = V 2/2 + Vi/3 - 0.248t for Vi/V 2 < 3/2, while the 
energy per site of the CDW II phase is -Ecdw ii = V 2 /2 + 
Vi/4. At t = 0, V\ = 0 both states have the same energy. 
At finite t, the CM phase has lower energy due to the first 
order quantum correction, which is absent for CDW II. 
Since Eqms grows faster with V\ than -Ecdwit there 
must be a crossing at a critical value of the interaction 
signaling a first order phase transition into the CDW II 
state. 


C. Cut C: Semimetal—CMs 

The cut at V\/t = 0.4, labeled cut C in Fig. 3j probes 
the transition between the semimetal phase and the CMs 
phase. In Fig. [6] we present the entanglement entropy S 
as a function of V 2 , the interaction that drives the phase 
transition. 

As explained in Sec. |III A| the entanglement entropy 
of the semimetal phase depends strongly on the bond 
dimension y. At V 2 /t ss 2.3 we observe a sharp transi¬ 
tion to a decaying entanglement entropy that does not 
depend strongly on y as V 2 is increased. However, the 
change in entanglement entropy is not as abrupt as in the 
CMs-CDW II case of the previous subsection as would 
be expected from a level crossing in a first order transi¬ 
tion. We can see a kink in S, but for decreasing V 2 it 
smoothly approaches the value of the semimetal phase. 
In addition, the energy is a smooth function V 2 which 
together with the previously mentioned behaviour of S 
suggests that the transition is of second order. The cor¬ 
relation length not depicted here shows a similar smooth 
behavior as in the transition between the semimetal and 
CDW I phases in Fig. [4] which provides further evidence 
in favor of a second order transition. 
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FIG. 6. Entanglement entropy at Vi/t = 0.4, labeled cut C 
in Fig. [3] probing the phase transition between the semimetal 
and CMs phases. Going towards the phase boundary from 
high V 2 , the entanglement entropy smoothly approaches the 
value in the semimetal indicating a second order transition. 


D. Cut D: CDW I Kekule-Semimetal-CDW II 


FIG. 7. Correlation length at Vi/t = 6, labeled cut D in 
Fig.0 In this plot, several phase transitions are made vis¬ 
ible. With increasing V 2 , the CDW I phase turns into the 
Kekule phase which itself goes into the semimetal phase. The 
transition from semimetal to CDW II can be identified by 
the peak in the correlation length. Inset: The vicinity of the 
Kekule-semimetal transition showing the discontinuity of the 
entanglement entropy as a function of V 2 


The next cut we focus on is labeled cut D in Fig. [3] at 
Vi/t = 6. The correlation length as a function of V 2 is 
shown in Fig. [7] Starting from low V 2 we first cross the 
phase boundary between the charge density wave CDW I 
and the bond-ordered Kekule phase. As expected for a 
gapped phase, the correlation length deep in the CDW I 
phase depends only weakly on the bond dimension. As 
we approach the phase transition, this behaviour changes 
since the many body gap is closing. The strongly in¬ 
creasing correlation length which peaks at V 2 /t « 1.5 
suggests a second order transition between the CDW I 
and Kekule phases. Moreover, we can see the presence of 
a remnant charge order on top of the bond modulation in 
the Kekule phase close to the critical point. This signals 
a slow onset of charge ordering as a finite size effect as 
we approach the phase boundary from the Kekule side 
which would be in contradiction to a level crossing for 
which the charge order should change suddenly. A sim¬ 
ilar second order transition was reported in the related 
antiferromagnetic Heisen berg model between a Neel and 
a plaquette phase jffl 43 l 44 l This transition was conjectured 
to show deconfined quantum critical behavior ! 45 ! 46 ! 

The fact that the correlation has not fully converged 
with x even in the center of the the Kekule region in¬ 
dicates that this phase is strongly entangled. This can 
be attributed to the fact that it is not a strong coupling 
phase in which the gap is determined by the interaction 
but rather by the hopping energy scale. As V 2 is increased 
we cross the boundary between the Kekule and the 
semimetal phase at V 2 /t w 1.9. From the data presented 
in Fig. [T]it is not possible to reach a conclusion about the 
nature of the transition. However, from the Landau the¬ 
ory perspective the semimetal-Kekule transition has to 
be first order for the following reason. The Kekule bond 


order is in the twofold representation E[ = 
and the lowest order non-trivial term in the Landau free 
energy is the scalar (F^) 3 — 3 E^iE'n) 2 which is of third 
order, therefore signaling a first order transition. A first 
order transition is also consistent with the entanglement 
entropy which displays a small discontinuity when cross¬ 
ing the phase boundary (see inset of Fig. [ 7 ]). 

In the semimetal phase, we observe an astonishingly 
small correlation length for a critical state. As explained 
in Sec. |III A| the hopping strengths around and along the 
cylinder are strongly renormalized. This shifts the Dirac 
nodes away form the K and K' points thus leading to an 
effectively gapped state for a finite cylinder circumference 
with small correlation length. 

Finally, at V 2 /t « 2.4 the semimetal turns into the 
CDW II. The diverging correlation length on the CDW II 
side of the transition points towards a second order phase 
transition. Moreover, the strong dependence of £ on the 
bond dimension suggests that the ground state is becom¬ 
ing critical as we approach the transition point, espe¬ 
cially compared to the behaviour of £ deep in the CDW II 
phase. On the other hand, from the semimetal side the 
correlation length is discontinuous, which could signal a 
first order transition. However, the latter observation has 
to be taken with care since the semimetal phase, being 
a critical state, cannot be represented faithfully with a 
finite bond dimension. 


E. Cut E: CDW I-Kekule-CDW III-CDW II 

The last cut we address is labeled cut E in Fig. [3] and 
fixes Vi/t = 9.2. The entanglement entropy is presented 
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FIG. 8. Entanglement entropy at Vi/t = 9.2, labeled cut E 
in Fig. [3] The smooth behavior of S at the phase bound¬ 
ary between CDW I and Kekule phase indicates a second or¬ 
der transition in accordance with the data of the correlation 
length in Fig. [7] At Vi/t ~ 2.8, the entanglement entropy is 
discontinous, which clearly indicates a first order transition 
between the CDW III and CDW II phases. 


in Fig. [8] With increasing V 2 we cross the second order 
phase transition from the CDW I to the Kekule phase 
at Vi/t ss 2.2 discussed above. The entanglement data 
shown in Fig.[8]supports this statement since S smoothly 
increases in the CDW I phase as we approach the phase 
boundary. 

At Vi/t ss 2.4 the system undergoes the transition to 
the CDW III phase. Unfortunately, a conclusive state¬ 
ment about the order of the critical point cannot be 
drawn from the present data; especially since the en¬ 
tanglement entropy in the CDW III phase still shows a 
significant dependence on the bond dimension. 

The last phase boundary in this cut lies between the 
CDW III and CDW II phases at V 2 /t « 2.8. Here, we 
can see a clear discontinuity of the entanglement entropy, 
similar to that in cut B (Fig. [ 5 ]). Moreover, we see a 
kink in the energy when crossing the phase boundary. 
Together these features indicate a first order phase tran¬ 
sition. 


V. DISCUSSION AND CONCLUSIONS 

The spontaneous emergence of the Haldane Chern in¬ 
sulator phas^ due to interactions in the half-filled spine¬ 
less honeycomb lattice model has bee n questioned by dif¬ 
ferent methods since its proposal! 9 19 12 — ' '— I In this work 
we have used the iDMRG method and have mapped out 
the phase diagram of this model. The algorithm allows 
for quantum fluctuations to play a role but minimizes 
the finite size effects due to its intrinsic infinite cylinder 
geometry. 

One of our main results is the absence of the Haldane 
Chern insulator state in this model. Even when initial¬ 


izing the iDMRG calculations with a time reversal sym¬ 
metry broken chiral wave function, the state did not re¬ 
main stable and evolved into the respective competing 
phases upon applying the algorithm. It seems therefore 
that quantum fluctuations indeed jeopardize emergence 
of the Cl phase and work against the naive mean field ex¬ 
pectations. Although there is no reason to believe that 
the situation is different for the 7r-flux model, which also 
hosts a mean field Chern insulating phase that is absent 
in exact diagonalizationf 19 Til there is more hope for other 
models with quadratic band point touchings. In particu¬ 
lar, the interacting kagome lattice hosts a Chern insulator 
state within mean field theory when the filling is tuned to 
the quadratic band point touching between the flat band 
and one of the two dispersive bandsJTH Its presen ce is 
predicted within the renormalization group approach 47 — 
which guarantees its robustness at sufficiently low inter¬ 
action strength. Adequate substrate engineering can also 
lead to such quadratic band point touching, potentially 
favoring the Chern insulator state 3^1 
In addition, we have theoretically accounted for two novel 
competing phases, the CDW II and CDW III. Their po¬ 
tentially large unit cells turns their identification within 
exact diagonalization or mean field studies challenging. 
On the one hand the CDW II phase stems from a de¬ 
generate strong coupling phase that has no first order 
correction in powers of t/V \Although the degener¬ 
ate manifold includes a distinguishable stripe-like phase, 
our iDMRG calculations find that a superposition state 
between classical ground states has lower energy. The 
particular superposition that is realized seems to be de¬ 
termined by the cylinder geometry where bond strength 
along and around the cylinder are inequivalent. The 
CDW III on the other hand, occurs in a finite region 
around the classical transition between the CDW I and 
II states defined by the line V\ = and has a clear 
twelve site unit cell. In this region, the first order cor¬ 
rections t/\ 1,2 lift the classical degeneracy and stabilize 
CDW III state. Both the semiclassical limit and our nu¬ 
merical iDMRG data show consistent charge structure 
factors confirming the semiclassical nature of the state. 
The CDW III phase presents as well a distinctive bond 
ordering that was understood from its charge order. 

From the evidence presented in this work we have 
established that the CMs state occurs from first order 
perturbation theory in powers of t/V 2 of an otherwise 
classically degenerate ground state. This part of our 
work establishes the existence and robustness of this 
single particle charge order beyond numerical approxi¬ 
mations. Our semiclassical analysis, corroborated by the 
iDMRG numerical evidence, provides an explanation of 
why this phase has a many-body gap of the order o f the 
hopping strength t/V 2 , as was previously noticed LL2120| 
rather than determined by V 2 . Both the semiclassical 
treatment and the numerical data we obtain establish 
that the CMs state has a finite subl attice imbalance, a 
feature that was still under dispute II21I2) Moreover, we 
have reported a characteristic bond order not addressed 
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in previous studies. As in the CDW II and CDW III 
cases such bond order is determined by the CMs charge 
order pattern. Our results show that the CMs is stable 
for small V\ and finite t. Increasing the former or 
reducing the latter leads to the CDW II ground state 
through a first order phase transition. 

Furthermore, we have established the existence of the 
Kekule bond order which occurs in a region of the phase 
diagram that appears to be smaller than that predicted 
by mean field and exact diagonalization. 

Finally, we have used the fact that the iDMRG method 
treats an infinite system as opposed to a finite cluster 
to analyse the different phase transitions using in 
particular the entanglement entropy S and correlation 
length £. The conclusions drawn from analyzing these 
quantities have been complemented by the continuous or 
discontinuous dependence of the ground state energy on 
the {V U V 2 } interactions. However, conclusions about 
the order of phase transitions from and to the semimetal 
phase should be taken with care due to the gapless 
nature of the semimetal ground state. 

To conclude we have established the phase diagram of 
spinless fermions hopping on the half-filled honeycomb 
lattice with nearest- and next-to-nearest neighbor inter¬ 
actions and characterized the phase transitions among 
the different phases with iDMRG. Our results provide 
solid evidence that Chern insulating phases are far more 
elusive than previously thought and so alternative routes 
are necessary to drive these kind of topological states 
from strong electronic correlations in general. 
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Appendix A: Strong coupling perturbation theory 

In this appendix we discuss the exact diagonalization 
method to determine the ground state in the strong cou¬ 
pling limit t -C Vi i2 . Consider splitting the Hamiltonian 
in Eq. Q into H = Hy + H t with 

H t = -tJ2( c l c j +h.c.) (Al) 

(*d> 

h v = v 1 J2 riirij + V 2 E "*"j- ( A2 > 

(ij) ((*J» 

In the strong coupling limit, we can obtain the ground 
state of H by diagonalizing Hy first and considering H t 
as a perturbation. The eigenstates of Hy 

Hy'lpn^m — (AS) 


where m accounts for degeneracies, are simple to compute 
because charge is conserved at every site in the absence of 
hopping. Hy is thus already diagonal in the occupation 
basis, i.e. 

I rn= n c !i°>’ ( A4 ) 

C 1 "’ Tn =l 

with = 0,1 the occupation coefficients for the 

n, m eigenstate. The classical ground state manifold is 
spanned by l^ 0 ’" 1 ), with m = 1,..., M. 

For small t/V-\ . 2 , we can disregard the classical states 
with n > 0, and project the H t Hamiltonian into the 
classical ground state manifold. Dropping the label n = 0 
from now on 

(Ht) mim2 = (i’ mi \H t \^). (A5) 

We can now diagonalize H t and select the eigenstates of 
lowest energy eo 

H t v a = e 0 v a , (A6) 

where a = 1 ,,D and D is the true ground state de¬ 
generacy. The final ground state of H is spanned by 

M 

\GS) a =5>? IVH ■ (A7) 

m= 1 

Obtaining the coefficients u™ is relatively simple because 
the effective size of the Hilbert space M is much smaller 
than the full size of the Hilbert space of H. To distinguish 
different phases, we recall that in finite size exact diag¬ 
onalization there is no spontaneous symmetry breaking, 
because all states related by symmetry are degenerate 
and will be present in the ground state manifold. In the 
simpler case when first order quantum corrections are 
zero, a phase can be characterized by inspection of the 
classical charge patterns of every eigenstate. With quan¬ 
tum corrections, however, the ground state can only be 
characterized by computing correlation functions evalu¬ 
ated in the ground state, from which order parameters 
can be obtained. Correlation functions for charge order 
parameters can be expressed in general as 

Pij... = tr (GS\ cjCjCjCj ... | GS ), (A8) 

which is given explicitly by 

fti... = £(C)XCrC7--- (A9) 

m,a 

The behavior of correlation functions can then be used 
to identify the different phases. We have used this 
method to diagonalize three different clusters: two 12 
site clusters with the periodicities of the CDW II and 
CDW III phases, and a cluster of 3 x 3 unit cells or 18 
sites, i.e. the L y = 6 version of the cluster in Fig. [l] The 
results can be summarized as follows. 
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In the case of the CDW II cluster, we find only two 
possible classical ground states with degeneracies 8 and 
2 for generic values of V 1 /V 2 . The projection of H t in 
these subspaces is zero for both cases, which means that 
both are stable strong coupling phases. The first can be 
identified with the classical version of the CDW II state, 
where the two values of the occupancies are 0 or 1. The 
8 states correspond to the 8 equivalent ways to ensemble 
the CDW II pattern in the given unit cell: two of them 
correspond to the stripe-like phase shown schematically 
in Fig. [3^c) while the other six correspond to those states 
with defect lines determined by rotated hexagons, as de¬ 
scribed in sec. |III| E. The second is the CDW I state, 
where the two states have a single sublattice (A or B) 
fully occupied. The energy per site of these states is 
given by 

-Ecdw ii = 2^2 + 4V1 p 4 < 4 (A 10 ) 

Ecdw 1 = ~y~ 4 < (AH) 

In the case of the CDW III cluster, for generic values of 
V 1 /V 2 we find two possible classical ground states with 
degeneracies 2 and 2. The first correspond to the stripe¬ 
like phase mentioned before, which is also commensu¬ 
rate with this cluster, and is therefore assigned the label 
CDW II. The second corresponds to CDW I. The ener¬ 
gies per site remain the same as those in Eqs. |A10|Alf| 
In addition, in this cluster there is another classical state 
with degeneracy 18 and energy E* = |T4 + gVi which 
requires consideration. This energy is always higher than 
Ectw ii or Ecdw i> except at the special point Vi = 414 
where the three cross, Eq dw ii = Bcdw i = E*. The 
reason why this state must be considered is that quan¬ 
tum corrections split it to first order in t , thus lowering 
its energy. Since neither CDW II or CDW I is affected by 
t to first order, there is a small region around Vj = 414 
where the energy of this state is lowest. Inspection of the 
phase diagram in Fig. [3] suggests that this is the CDW III 
phase. The ground state of this phase has degeneracy one 
and energy 

5Vj Vj 

Ecdw hi = -g- + -g— 0.236f p 4 ~ 4 (A12) 

To confirm the nature of this state, we have computed 
the Fourier transform of the correlation function C(q) = 
J2ij e^ Xi ~ x ^ q pij, and confirmed that it compares favor¬ 
ably with the Fourier transform of the charge pattern of 
CDW III shown in Fig. |gd). 

In the case of the 3x3 cluster, at t = 0 we find four 
possible classical ground states as a function of increas¬ 
ing Vi, with degeneracies 234, 108, 36, 2. The projection 
of H t into these classical ground states is finite only for 
the first two, and the pattern of degeneracies at low en¬ 
ergies becomes (4,10,4,...) for both. This low energy 


pattern is the same as the one that is obtained for the 
CMs phase in exact diagonalizationj 4 ^ pointing to the 
fact that these two states represent the CMs phase (fur¬ 
ther evidence to support this claim is also shown below.). 
The third state is an intermediate state with the same 
energy as the CDW III classical states, but which is not 
corrected by quantum fluctuations at any Vj/Vjj. The 
fourth is again the CDW I. The energies per site of the 
three phases in the 3x3 cluster are given by 

( Vk 1 Vi _ q 940^ Vi < 3 

£ ™- = { A + A-0.055* |<!i<3 (A13) 

5 Vo Vj 

E* = -^ + -± 3<^<4(A14) 

3Vn 

Ecdw 1 = ~2~ 4 < p| (A15) 

As explained in the main text, the energy per site of 
the CDW II phase Ecdw ii is lower than the classical 
energies of any state of the 3x3 cluster. Since the 
intermediate state does not have quantum corrections, 
its energy satisfies Ei nt > E cdw ii and it is never re¬ 
alized. The CMs state, however, has Ecms < -Ecdw ii 
for small enough Vj once the quantum corrections are in¬ 
cluded. Further evidence that this is the CMs state can 
be obtained by computing correlation functions. Fourier 
transforming the real space indices ij ... and taking the 
combinations defined in the main text for the representa¬ 
tions B 2 and G, we can compute any correlation function 
of order parameters {GS\ O 1 O 2 ■ ■ ■ IGS 1 ) with O = G, i? 2 . 
The simplest scalars that signal the presence of charge 
order are simply the two point functions 

Sj = Bf, (A16) 

5 2 = |G 1 | 2 + |G 2 | 2 . (A17) 

We have computed these correlation functions in the 
ground states of the 3x3 cluster with V\ < 314 and 
finite t and found that they are both finite. However, 
this is not enough information to distinguish the CMs 
phase. In order to detect its particular order we have 
also computed the correlation functions of higher order 


for this ground state 


S 3 = Im(Gl - 3G 2 G 2 ), 

(A18) 

S 4 = B 2 Im(GiG 2 ), 

(A19) 

S 5 = [Im(G 1 G' 2 )] 2 , 

(A20) 

S 6 = B 2 Re(Gl - 3GiG 2 ), 

(A21) 


and found that only Sq is finite, therefore confirming the 
presence of the CMs phase with the strong coupling ap¬ 
proach for Vj < 314. As for the CDW III we have also 
compared the semicalssical and numerical structure fac¬ 
tors and found good agreement, supporting the semiclas- 
sical interpretation of the CMs phase. 
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